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I. INTRODUCTION 



The recent progresses of experimental studies of Bose- 
Einstein condensates (BECs) in optical lattices (OLs) JTj 
and in particular, the direct observation of a gap soli- 
ton have stimulated many studies of the dynamics 
of matter waves in periodic media (see e.g. the reviews 
and references therein). In relation with nonlinear 
matter waves, OLs provide a tool for changing the effec- 
tive properties of the atomic medium allowing the exis- 
tence of solitary waves in condensates with positive scat- 
tering length. Moreover, deformed OLs offer additional 
possibilities for manipulating the soliton dynamics. In 
particular, smoothly modulated lattices are effective in 
controlling the quasi-one-dimensional dynamics of small 
amplitude gap solitons [1,0. 

The study of various types of localized defects and their 
interaction with matter waves has been considered in 
different contexts, ranging from experimental studies Q 
and different theoretical works studying the emission of 
vortices in the superfluid flow of a BEG around a de- 
fect 1^, the generation of gray solitons and sound waves 
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by a moving defect Q , the drag force appearing at large 
velocities of a moving defect in a BEG (81, and the com- 
bined effect of a static defect and trapping potential . 

In Ref. we have explored theoretically some fea- 
tures of the stationary configurations of a BEG under 
the simultaneous effect of an OL and of a strongly lo- 
calized defect. This leads to the existence of localized 
states, termed nonlinear defect modes, which can be 
driven through hundreds of periods of an OL. In this 
paper we extend the analysis of Ref. addressing 
the problems of existence, stability, and generation of 
stationary modes localized in the vicinity of the defect, 
and consider the possibility of the dynamical interac- 
tions between defect modes. We provide details of an 
approximate analytical description of defect modes, out- 
lined in [lol | and report thorough numerical simulations 
of their dynamics in order to understand their stability 
properties. To do so, we will exploit the fact that in the 
presence of an OL, the wave function of the BEG can be 
expanded on a set of Wannier functions (WFs) 0, 0, 
localized on each potential well and modulated by an 
envelope. 

Before going into details, we would like to mention that 
the theory of nonlinear defect modes in one-dimensional 
(ID) periodic media, is an issue of high relevance also 
for the theory of photonic crystals. The case of lattices 
with shallow grating was recently addressed in p3j within 
the framework of the coupled-mode equations for the 
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counter-propagating waves. In the same approximation, 
hght propagation in nonhnear fibers with slowly modu- 
lated Bragg gratings was considered in Ref. jl^- The 
research reported in the present paper, as well as in ear- 
lier publications ^ .10] , differs not only by the physical 
model, but also by the fact that the grating is consid- 
ered to be deep and the dynamics is studied numerically 
within the full ID equation. 

In this context it is also relevant to mention studies of 
impurity modes in the framework of the discrete nonlin- 
ear Schrodinger equation |l6l |. since in the tight-binding 
approximation, i.e. in the limit opposite to shallow lat- 
tices, the nonlinear Schrodinger (NLS) equation with a 
periodic potential can be mapped into the discrete NLS 
equation |l[ll[l2. 

The structure of this paper is as follows: In Sec. Ullwe 
expand the wave function in terms of the WF's and derive 
coupled equations governing the evolution of the enve- 
lope. Then, in Sec. IIIII wc provide the simplest analytical 
solutions corresponding to stationary defect modes. We 
also study the dynamics and stability of these modes by 
direct numerical simulations of the ID Gross-Pitaevskii 
(GP) equation. Another interesting question is the pos- 
sibility of moving these defect modes, which we consider 
in Sec. II VI In Sec^we propose several ideas for the gen- 
eration of defect modes in realistic experiments. Finally, 
in Sec. IVIjwe summarize the outcomes. 



II. 



STATEMENT OF THE PROBLEM AND 
ANALYTICAL CONSIDERATION 



A. Bloch and Wannier function expansions 

We consider a BEG trapped in a quasi-lD configura- 
tion, i.e. tightly confined by a trap potential in the trans- 
verse direction, providing a cigar-shaped geometry. We 
also suppose that the system is subject to a ID periodic 
potential, which is created by laser beams and whose axis 
is directed along the longitudinal direction of the trap. 
Provided the density of atoms is small, such a system is 
governed by the ID GP equation (see e.g. 01 [l3 for the 
details of derivation), which in dimensionless variables 
reads 



.dtp d'^ip 



dt 



+ Vi{x)4' + Vd{x)4' + (j\^\''4>. (1) 



Here a — \ stands for repulsive and a = — 1 for attrac- 
tive atom-atom interactions, Vi{x) is an OL potential and 
Vd{x) is a defect potential. The dimensionless variables 
are chosen so that the energy is measured in the units 
of the recoil energy Er = h^k'^ / {2nn) where k = ii/d and 
d is the lattice constant. Then the temporal and spa- 
tial coordinates are measured in the units of E^/h and of 
d/TT, respectively, and the period of an OL potential is tt: 
Vi{x -\- tt) = Vi {x) . The dimensionless macroscopic wave 
function is normalized as follows J jV'prfa; = ^TiNk\as\, 
where is the scattering length. 



In order to provide a formal derivation of the approxi- 
mate equations we will assume that the defect potential 
is localized in space on a distance smaller than one lat- 
tice period, what mathematically can be expressed as 
I T^^E^ 1*^ 0{'K^'^) (the numerical studies presented bel- 
low show that this region can be significantly enlarged). 
For the sake of simplicity, we will restrict our analysis to 
symmetric OL potentials Vi{x) = Vi{—x) and represent 
the defect potential in the form Vd{x) — V+{x) + V^{x), 
where V±{x) = ±V±(— x) are symmetric, V+, and anti- 
symmetric, V-, parts, respectively. 

In the absence of nonlinearity, Eq. ^ describes the 
well-known problem of a quantum particle in a periodic 
potential which arises, for instance, in description of the 
interaction of an electron with a defect in a ID crystal 
lattice. We will make use of this fact and use approxi- 
mation techniques which are well developed in solid state 
physics 13 • 

To this end we introduce an expansion of the macro- 
scopic wave function 



oo „i 

1p{x,t) = '^ dk(pak{x)fa{k,t) 
^— n — 1 



(2) 



over the orthonormal set of the Bloch functions ipakix) 
of the linear eigenvalue problem 



d^ 



dx^ 



+ Vl{x)(pak = Ea{k)ipak- 



(3) 



Hereafter a stands for the number of the zone, subject to 
convention that the lowest allowed band has zero number, 
k is the wavevector considered in the first Brillouin zone 
(BZ), k G [—1,1]. Then it is a straightforward algebra to 
obtain from Eqs. Q-© the following equation for the 
envelope faik, t): 



iU{k,t)~ E^{k)fo,{k,t) 



dx^pakix)\^|J\'^■lp{x,t) 



+ J dk' J dx(pak{x)Vd{x)ipa'k'{x)fa'{k' ,t). (4) 

a' 

In this paper we drop the limits in all integrals, taking the 
convention that integrals with respect to x are considered 
over the whole line and integrals with respect to k are 
considered over the first BZ, i.e. / dx... = dx... and 

J dk... = J^^ dk.... 

Now we introduce WFs 



Wo 



^ix) 



1 

71 



dk (pak{x)e' 



(5) 



which make up a complete orthonormal basis and can 
be chosen real and exponentially decaying Each 
WF is centered at x — nn, where n is an integer. The 
Bloch functions can be obtained from the WFs from the 
inversion of Eq. (|SJ) through the relation 
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^{x)e' 



(6) 
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Introducing the notation 



(7) 



and defining the Fourier transform of V^^\ : 

V'^'kk' \ T /-nn'inik' n' —kn) /o\ 



which possesses the evident symmetry V^l^, = Vj^/^, we 
rewrite the defect term in Eq. Q in the form 



In Eqs. ijHl and (j^ we have used the convention that 
ah sums with respect to n and a are considered over 
ah integers and all nonnegative integers, respectively, i.e. 

n ■■■ = 2^n=-oo ■■■ •■■ = 2^Q=0 ■■■■ 



J J dx(pak{x)Vdix)ipa'k'ix)fa'ik' ,t) 

a.' 

= E/ dk'V^^,U'ik',t). (9) 



Next we rewrite the nonlinear term in Eq. as 



E 



^ J dk.dk^dk^wziZl 
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71711712^3 i7T{k2'n2'\-k^n^—kn—kini) 



fa^ikl,t)fa^{k2,t)fa^{k3,t) 



ni7i2n^ 01O2O3 



r 



where 



dxWart^aini'^a2n2 ^03 



(11) 



In writing Eq. (|10|l we have used the property that 

W^Zi^^laf = W^^^^cT.Sf^^""'"'"" (see P 113 for more 
details) and have computed the sum over n. Now that 
we have rewritten Eq. ^ in an appropriate form for 
our analysis we will discuss the approximate analytical 
description to be constructed in this paper. 



B. Constrains and approximations 

We will be interested in matter-waves of small ampli- 
tude, which are well described in the effective mass ap- 
proximation [^Il9l|. and thus having spatial extension. A, 
much bigger than the lattice constant: A ^ tt. Then the 
respective characteristic scale Ak ^ 2tt/X of fa{k,t) in 
the momentum space is much smaller than the vector of 
the reciprocal lattice (since Ak • A ^ 1) 



Afc < ksz = 1 



(12) 



Meantime the scale of the variation of V^j^, with respect 
to either of the arguments exceeds one. In other words, 
we assume that fa{k,t) is a strongly localized function, 
compared with Vj^^, , the latter being weakly dependent 



on both k and k' . Assuming also that /^(fc, t) is localized 
about some point feg in the first BZ (i.e. |A:o| < ksz) we 
can approximate 



v5lo.'{k\t)dk' ^V^'J I U'{k',t)dk' 



U>{k',t)dk'. (13) 



WFs are known to be exponentially localized in 
space 01 and the "radius" of localization decreases as 
the potential depth is increased, assuming the period re- 
mains constant (for examples see e.g. jsl ll^ ) . 

Subject to these conditions and taking into account 
the inequality p2ll and the fact that f{k,t) is localized 
about fco we conclude that the exponential factors e^^i^i'^ 
in (|10|l can be substituted by e™^^°'^ (with exponential 
accuracy), what allows us to introduce 



l"2n3 pi(n2+n3—ni)ko-K 
ia2a3 



(14) 



nin2n3 

and rewrite Eq. p(J|l in the form 

E W^aaia2a3(7Qi(fcl>^)/a2(^2,i)/a3(fc3,i)), (15) 

010203 

where, for the sake of brevity, we have introduced the 
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notation 



{a{h)b{k2)c{h)) = 




dkidk2dk'i5k+kiM+ki 



xa{ki)b{k2)c{k3). (16) 

Whenever V^j^? ^ the defect term (O resuhs in 
tunnehng of atoms between the bands a and a'. The 
nonhnearity (|14() is another cause for inter-band tunnel- 
ing. Therefore one can expect that a one-band approx- 
imation, successfully reproducing gap solitons in some 
homogeneous lattices [ll|, may not work properly when 
the combined effect of the defect and the nonhnearity are 
considered. Below we will see, however, that this is not 
always so. 

We will also require that the potential Vi[x) is such 
that the lowest even (including zero) and odd bands to 
have even and odd WFs, what, in particular, means that 

Wan{x) = {-l)°'Wan{-x) . (17) 

This is the case of a symmetric OL potential, Vi{x) = 
Vi{—x) having a local minimum at x — 0: Vi{x) > V;(0). 



The imposed constrain is not strong, and in particular is 
satisfied by the cos-like potential [see (|?7| below] model- 
ing most of the presently available experimental settings. 

Under the imposed conditions, one immediately ob- 
tains from (|nj that W^™^^^.^ = whenever a -I- ai -|- 
Q!2 + ct3 — 2p-|-l withp being an integer. This means that 
the lattice originates an effective coupling of the bands 
with the same parity (i.e. for example atom transfer be- 
tween band a = and a = 2 or between bands a = 1 and 
a = 3). This fact leads us to the next simplification of 
the problem, considering the two-band approximation as- 
suming that initially all atoms possess energies in the two 
lowest bands, and thus only the bands (with a — and 
a = 1) will contribute to the BEC dynamics. Thus we 
will impose that fa{k, t) for a = 0, 1 and fa{k, t) = 
for a > 2. Since this assumption is not supported by any 
rigorous statement (and the defect itself can result in the 
coupling of the bands) below the respective results will 
be compared with direct numerical simulations. 

Then, taking only a, /? = 0, 1, /3 7^ a we obtain the 
final expressions for the nonlinear term 



W^a{fo,{kl,t)fa{k2,t)Uk3,t))+W^P (2 (/^ (fcl , i)/„ (fcz , i)//3 (A^S , t)) + (/a (^1 , i)//3 (^2 , t)//5 (^3 , *))) , (18) 

Now defining 14 = 27rV^^'=«, V = 2'KV^f° = 2ttV^°''° and M4/3 = with Wom = 1^0001 = (the last two 

approximations are made on the basis of strong localization of the WFs resulting in the estimate |M^aQ^Q2a|| <C 
l^aafaaaal least One of Uj ^0, See the examples in Table below) we obtain the equations 

ifo(k,t)-Eo{k)fo{k,t)-^ J fo{k',t)dk' J fi{k\t)dk' ~ aWoo{7oiki,t)Mk2,t)Mk3,t)) 

-aWoi [2(7i(fci,t)/o(A:2,t)/i(A;3,t)) + (7o(fci,i)/i(fc2,i)/i(fc3,t))] =0, (19a) 

i/i(fc,t) - E^{k)h{k,t) - ^ J h{k',t)dk' - ^ J fo{k\t)dk' - aWn{7iiki,t)fi{k2,t)f\{k3,t)) 

-aWio [2(7o(fci, i)/i(A;2, t)/o(A:3, <)) + (7i(fci, i)/o(fc2, t)/o(fc3, <))] = . (19b) 



In the vicinity of the maxima of the functions fa{k, t), 
i.e. in the vicinity of fc = fco, one can expand 



Ea{k) = Ea + Va{k - fco) 



{k - fco)' 



(20) 



where Ea — Eaiko), Va 



_ dE^(k) 
~ dk 



k—ko 



2Ma 

is the group ve- 



locity of the mode a, and Ma — y^-^§T^ 
effective mass. 

C. Dynamical equations 



is its 



k—k(i 



Multiplying Eqs. (|19|) by exp[i(fc — fco)a;] and integrating 
it with respect to k over the BZ, we obtain a system 
of coupled nonlinear Schrodinger (CNLS) equations with 
delta impurities 



Let us now introduce the function 

Ux,t)^ [ e*('=-'=»)V„(fc,t)dfc. 



(21) 
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1 



dfo . dfo 

* dt dx ^ 2Mo dx^ 



' dt ^ 9a; ^ 2Mi 5^2 



- So/o - dWool/oP/o - aW^oi (2I/1IV0 + /i /o) = Vo'5(x)/o + V5{x)h, (22a) 

- - <7Wii\h?h - <yW„i (2I/0P/1 + /o/i) = T>'5(.t)/o + (22b) 



These equations must be consistent with the expansion 
used for their derivation which imphes that the effective 
masses should be of order of one, \Ma\ ~ 1. These condi- 
tions are satisfied in reasonably large ranges of parame- 
ters for the most interesting potentials from the physical 
point of view, when the depth of an OL is of the order of 
a few recoil energies. 

To recover the wavefunction ip{x,t) from the solutions 
of Eqs. we make use of Eqs. lO and ^ together 

with the definition pifl and obtain 



ip(x,t) 



V2 



e'''''"'w^nix)U{mr,t). (23) 



In what follows we concentrate on the modes bordering 
the gap, i.e. modes with fco = ksz and with Va = 0. In 
that particular case and taking into account that we are 
working in the two-band approximation, Eq. H23f) takes 
the simple form 



V2 



^wori{x)fQ{m:,t) +win{x)fi{mr,t) .(24) 



Now we also have that, Va and V are originated by the 
symmetric and antisymmetric components of the defect. 



I.e. 



Va = ^Y.^-IT+''' j V+{x)Wan{x)Wan-{x)dxi^^&) 
^ = ^E(~l)"^"' f V-{x)won{x)win'ix)dx{25h) 



Indeed, using that Wan{x) = Wa{x — 7rn), i.e. the fact 
that WF's depend on (a; — 7rn), as well as the symmetry 
of WFs (113, we get 



V{x)[Wan{x)Wa'n'{x) + Wa-n{x)Wa' -n' {x)]dx 

[V{x) + {-ir+'''V{-x)]Wan{x)Wa'n'{x)dx (26) 



from which Eas. l25|l follow. 

Eqs. imply that even defects couple only modes 

of the same parity (and in the two-band approximation 
does not result in any coupling), while the odd defect 
couples bands with different parity (i.e. the two adjacent 
bands in the two-band approximation). 



D. Lattice potential and numerical test of the 
approximations 

During the derivation of Eqs. (|22l) a number of ap- 
proximations have been used and several constrains have 
been imposed. That is why before going into details of 
solutions we will test the reliability of the analytical ap- 
proximation. For the numerical calculations we will use 
the periodic potential in the standard form 



Vi{x) = Acos(2a;) 



(27) 



where A is an amplitude of the lattice potential (we recall 
that it is measured in the units of the recoil energy). 
By varying the amplitude of the OL one can control the 
magnitude of the parameters in Eqs. (|22|l . 

The first of our approximations has been to assume 
strong localization of the WF's. We have also used an ef- 
fective mass approximation which describes the behavior 
of the wave function near the boundary of the BZ what 
corresponds to the case of a small amplitude wavepacket. 
Thus, our theory is expected to be applicable for small 
amplitude wavepackets in relatively large amplitude OLs. 

We have calculated the coefficients in Eqs. (|22|) for 
different amplitudes of the OL, which are presented in 
Table m 25]. As one can see from the table, the non- 
linear coefficients in Eqs. (|22f) are of different orders of 
magnitude. Specifically, Wa/i ^ Waa ce ^ f3, what is 
a natural consequence of the symmetry of the WFs. In 
fact, this is a property of a large class of potentials rather 
than a feature of our specific model (|27|l . This approxi- 
mation allows us to simplify significantly the problem by 
neglecting the cross-phase-modulation terms, i.e. terms 
with 11^01 in Eqs. 

^From Table U we also see that in most of the cases 
the effective mass is indeed of order of one, as it was 
supposed in Sec. El The explicit forms of the WFs for 
the cases considered in the table can be found in Ref . 0| , 
where it is shown numerically that the WFs are indeed 
well localized on several (in the case A = —1; —3) periods 
of the OL, or even on a single period when A = —5. 



III. STATIONARY IMPURITY MODES 



A. One-band approximation. Bright and dark 
defect modes 



An important effect introduced by a localized poten- 
tial is the coupling between different bands. When the 
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TABLE I: Effective mass, energy, and overlapping coefficients at the boundary of the BZ and at the edges of the first gap for 
different amplitudes of the OL. In parenthesis we indicate the signs and numbers of identical hopping integrals contributing to 
the nonlinearity H^Qc«ic«2a3 defined in 11411 . 



Lowest band, q = 



A 


Mo 


Eo 


W^oo 


AAAA 

TT/UOOO 

'^''0000 


(1+) <ooo(6+) 


<ooo(8-) 


<ooo(8+) 


W^oooo(6+) <o°o'o'(12-) 


-1 
-3 
-5 


-0.163 0.471 
-0.85 -0.733 
-2.5 -2.076 


0.2505 
0.2983 
0.333 


0.3753 
0.5479 
0.6479 


0.0082 
0.0006 
0.0001 


-0.0064 
-0.0052 
-0.0022 


0.0035 
0.0005 
0.0001 


-0.0026 
-0.0001 
0.0000 


-0.0011 
-0.0000 
-0.0000 














First band. 


a = 1 








A 


Ml 


El 






W«>ff(l- 




W^iT/(8-) 


M/iTi'(8+) 


W^i'//f(6+) W^i%^(12-) 


-1 
-3 
-5 


0.1 

0.21 

0.3 


1.467 
2.166 
2.5 


0.2069 
0.2081 
0.2016 




0.2404 
0.3236 
0.3976 


0.0358 
0.0211 
0.0075 


-0.0163 
-0.0036 
0.0048 


-0.0008 

0.0035 

0.0025 


-0.0148 
-0.0079 
-0.0020 


0.0064 
0.0035 
0.0008 


Cross terms 


A 


W^oi 


11/0000 


T,r0011 

V^OOll 


Tj/OOOl 


Ti/0002 w0012 
I'^OOll ''''0011 


TI/-0022 M/OlOO 
''fOOll '''^0011 


TI/OlOl 

''^0011 


TI/0102 
''Vooii 


TI/0112 M/0122 
''''OOll '''^0011 


-1 
-3 
-5 


0.1306 
0.1091 
0.1224 


0.1815 
0.2221 
0.2703 


0.0399 
0.0105 
0.0028 


-0.0031 

0.0089 

0.0100 


0.0027 0.0073 
0.0049 0.0026 
0.0027 0.0005 


0.0026 0.0074 
0.0007 0.0019 
0.0001 0.0004 


-0.0063 
-0.0018 
-0.0004 


-0.0017 
-0.0003 
0.0000 


0.0017 -0.0029 
0.0004 -0.0001 
0.0001 -0.0000 



defect has a well-defined parity, the problem is further 
simplified, as it follows from (|25|l . Concentrating, first, 
on a defect of even symmetry, where Vd(2;) = V+(a;) we 
obtain that F = and the system of equations H22(l be- 
comes decoupled transforming in (a. = 0, 1): 



' dt 2M„ 9x2 



(28) 



This means that for even defects we can assume validity 
of the one-band approximation which was the approxima- 
tion used in Ref. [lU^ . 

The problem of interaction of a soliton of the NLS 
equation with a localized impurity has been thoroughly 
studied in literature (see e.g. Ref.[20j and references 
therein). Those results however are not directly appli- 
cable to our case, since fa is only an envelope of the 
WF's while the solution of our problem Q is given by 
Eq. (113) [or Eq. 

Since the first-band WFs are odd with respect to 
the spatial variable and the defect potential is localized 
around the zero of the WFs, for the case of the strong lo- 
calization |Vi| is negligibly small compared to |Vo|. Then 
the defect has much stronger influence on the modes of 
the zero band than on those of the first band. Therefore 
the consideration in this section is restricted to the low- 
est band, a = 0. Bright localized modes in this case were 
considered in Ref. [l3- Aiming to extend those studies 
by including new dynamical properties of the modes in 
the present subsection we recall the relevant formulas, as 
well as extend the analysis by including dark modes in 
the consideration. 

We will look for stationary modes in the one-band ap- 



proximation of the form: foix^t) 



-iEt 



00 (x) where 



A description and classification of bright defect modes of 
the NLS equation can be found in Ref.|23]. Below we 
recover some of the known results, which in our case are 
applicable to the envelope. Besides we will also discuss 
dark defect modes, as well as the effect of the negative 
mass and of the fine structure, i.e. WF's constituting a 
background of a defect modes, on the mode dynamics. 

A relevant parameter of the problem is a detuning of 
the mode energy outwards the band edge, which is de- 
fined SiS eo = E — Eq. Then considering Eq <^ Ex — Eq, 
what means smallness of the detuning, and keeping only 
the terms of the leading order, we can simplify Eqs. I|22() 
reducing them to the stationary NLS equation with a 
delta impurity 



1 <P 



1^0 



2Mn dx^ 



£000 - crW^oo^o = Vo(5(a;)0o- 



(29) 



Stationary bright modes of Eq. H29I) can be excited 
with the energy belonging to the gap, i.e. for a positive 
detuning: eo > 0. Introducing cry =sign(Vo) and taking 
into account that Mo < 0, one finds that the cosh-mode 
can be excited only if cr > and in a general case has the 
form 



00 (a;) 



^/2WWc 



00 



cosh [V-2Moeo(|a;| + xq)] 



(30) 



where 



4>o{x) is a real function and E is the energy of the mode. 



xo = ^ =atanh./— , 

with eo > e* and e* — — ^MoVq ■ the case Vq > 
the cosh-mode has one maximum and otherwise it has a 
two-hump profile (see Fig. 2 and the respective discussion 
in ReL 10]). 
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In the case of attractive interactions, i.e. when cr < 0, 
and only if V^q > one can find another bright locaHzed 
solution of Eq. - the sinh-mode 



00 



sinh [V-2Moeo(|a;| +a;o)] 



where 



1 



\/-2Afo£o 



atanhi / — . 



(31) 



(32) 



and Eq < e» [for examples see below Fig.n(a)-(c)]. 

Finally, in case of attractive interactions, cr < 0, con- 
sidering energy detuning towards the allowed band, i.e. 
Eo < 0, one can find a dark localized mode: 



£0 



tanh 



where 



cry 



00 



latanh 



v/£oMo(|x| + Xo) 








(33) 



(34) 



It follows from the obtained formula that for Vb < one 
can construct dark modes with an envelope having one 
hole at the defect position [see Fig.n(d)-(f)] while in the 
case Vo > the envelope of the dark defect mode has two 
holes [see Fig.n(g)-(i)]. 



B. Numerical results: Comparison with the 
analytic approximation and stability 

In what follows we will compare the previous analytical 
results with direct numerical calculations of the defect 
mode profiles. To do so we take the same form of the 
defect potential as in 



27r£ 



(35) 



This defect is characterized by its amplitude 77, width £ 
and the location of its center Xd- Unless otherwise speci- 
fied we will assume that Xd coincides with the minimum 
of the potential located at x = (this also implies that 
A < 0). The chosen defect potential can be generated 
by illuminating the condensate with a localized beam 
with gaussian profile and orientation perpendicular to 
the beams generating the OL potential. It is also easy to 
obtain from Eq. (|35|l the limit of a point defect (i.e. the 
delta potential) by taking £ — » 0. 



As it follows from the list of the hopping integrals W 
presented in Tabled for obtaining Woo with accuracy of 
about 3%, even in the case of amplitude ^ = — 1 it is 
enough to take into account only nonlinear overlapping 
integrals between neighbouring sites [see the definition 

The numerical study of cosh-modes was done in our 
previous paper ^3 ■ Here we complete analysis of bright 
defect modes by presenting in Fig. n(a)-(c) sinh-modes 
calculated analytically from Eqs. (|24|) and H31|) and com- 
pare them with modes calculated as numerical solutions 
of Eq. for different defect parameters. Considering 
the lowest band, ci; = 0, we find that sinh-modes exist 
only when > and eo < £*• In this case, the defect 
should be taken sufhciently large to ensure that differs 
appreciably from zero. 



In Fig. ^ (d)-(i) we present dark- modes in the form 
of tanh- function calculated analytically from Eqs. H24|) 
and (|33|) . Dark modes have energy shift Eq towards the 
allowed zone and exist at the boundary of the BZ only 
for attractive nonlinearity, cr < 0. We also observe that 
dark modes exist for any ratio Eo/^*- For negative Vd 
one can construct one-hole tanh-mode [see Fig.n(d)-(f)] 
and for positive Vd tanh-modes with two holes [Fig. 
(g)-(i)]. The holes however refer to the envelope of the 
WFs and therefore are not always visible in the profile, 
as it happens for instance, in Fig.n(h), (i). 



We have checked using analytical formulas as well as 
direct numerical simulations that by decreasing the de- 
fect width by factor of ten, the shapes of defect modes do 
not change appreciably and their amplitudes differ only 
by 1%, what confirms the good accuracy obtained from 
the delta-approximation for all £ < tt. Moreover, Fig. 
shows, that in the outlined parameter region, our ana- 
lytical approach gives the mode profiles with sufficiently 
high accuracy, even in the one-band approximation. Ob- 
viously, increasing the value of the defect strength, 77, 
the analytical approximation becomes worse. When 77 
becomes of order one, the difference between the analyt- 
ical predictions and the numerical results becomes ap- 



preciable. This occurs because such a defect originates 
modes localized on very few lattice periods and the ef- 
fective mass approach [ly fails. Then one expects that 
the tight-binding approximation [TH Il2|| would be more 
appropriate [2^ . 

The stability of the defect modes was tested by com- 
puting the evolution of the cosh-, sinh- and tanh-modes, 
starting from their analytical expressions and shown in 
Fig. [21(b), (d), (f), (h) and (j). Since the respective data 
are not exact solutions of Eq. (^), such initial conditions 
automatically introduce perturbations of order of a few 
percents with respect to the true stationary mode. In all 
the cases depicted in Fig. [3 one observes stable evolution 
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FIG. 2: In (a), (c), (e), (g) and (i) the initial (thick, black) and final (thin, red) profiles of the defect modes and in (b), (d), 

(f) , (h) and (j) their corresponding dynamics are shown. The initial profiles for cosh- modes is calculated by using the following 
parameters: A = —1, a = 1, £ — 0.1 and i] — 0.1, eo = 0.03 for (a) and rj — —0.5, eo ~ 0.09 for (c). The initial profiles in (e), 

(g) and (i) are the same as the analytical solutions depicted in Fig. 0(b) (f), (g), correspondingly. 




FIG. 1: Stationary defect modes of the lowest band {a = 0) for A — —1 and £ = 0.1. Thick and thin lines correspond to 
numerical and analytical solutions respectively. In (a)-(c) we show a sinh-mode for a = —1, Vd > 0, in (d)-(f) we plot the 
tanh-mode for a = —1, Vd < and in (g)-(i) two-hole tanh-mode for a = —1, Vd > are presented. In all the figures A, £ and 
r\ are given parameters of the lattice and of defect, while e,, Vb, and are computed through their definitions given in the 
text. The dashed lines show combined potentials whose scale is indicated on the right axes. 



of the modes. 



C. Coupled modes for an odd defect 



Coupling of the bands can occur either due to the 
nonhnearity or due to antisymmetric part of the defect, 
V- . In the present subsection we concentrate on the last 
possibihty, while still neglecting the cross-phase modula- 
tion. More specifically, we consider stationary solutions 
fa{x,t) = e~'^^^ (paix) , for which the system H22|) is re- 



duced to 



1 d^cbo 



+ £000 - crWQo(l)l = VS{x)(f>i , (36a) 
+ ei(t)i - (jWii<l>l = VS{x)4>o . (36b) 



For obtaining a stationary localized nonlinear defect 
mode of (|36|l . one has to take into account that in the 
case at hand MqMi < 0. Then solving the system (|36|) 
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one can compute 



0/3 (a 



cosh [6q(|x| — , 



sinh[6^(|a;| - xp)] 



(37a) 



(37b) 



where — y/2\Mje^\, the indexes {a, (3) are (0,1) for 
CT > or (1, 0) for (7 < 0, correspondingly, and the con- 
strain MaSp > is verified. The parameters x^^p are 
found from the following system 



tanh(6Qa;Q 



sinh{bpxp) 
cosh(6Q,a::Q,) 

cosh(6Q,XQ) 



V spM^Wo, 



2W, 



130 



tanh(6^a;;5) sinh(5^a;^) 



V 

e/3 



SaMfiWpp 



,(38a) 



(38b) 



Strictly speaking, mathematical consistency of the ob- 
tained system of equations requires smallness of both de- 
tuning parameters Eq and ep. This follows from the im- 
posed condition of the smoothness of the envelopes (jia.p 
compared to the lattice constant (what allowed us to in- 
troduce the effective masses). That is why, the applica- 
bility of the two-level approximation is determined by the 
width of the gap, A = IeqI + which should be small, 
in order that both detunings \Sa,i3\ be small. This can 
be provided by small amplitude OLs. On the other hand 
one should have WFs localized on a few periods of the 
OL what can be satisfied by increasing the amplitude of 
the OL. These arguments make applications of the sys- 
tem H37f) to be rather restricted and requires a numerical 
verification of the obtained solutions. 

To this end we choose the odd defect in the following 
form 



Vd{x)^V-{x) 



2Tri 



572 



(39) 



which according to Ea. (|25bp gives only antisymmetric 
component V of the defect while symmetric components 
Vo and Vi are zero. Next we recall the results presented 
in Table ^ where the effective masses for the potentials 
A = —1 were found to be of order of 0.1 while the width 
of the lowest gap is approximately 1. Thus, one can ex- 
pect that in this case the above requirement are approxi- 
mately satisfied. In Fig. Owe show an example of a defect 
mode, which first was computed using the analytical ap- 
proximation from the expressions H37a() . 1)3 7bl) and H24(l 
(i.e. the expansion over the WFs [see Fig. 13 (a)] in the 
potential shown in the panel (b). Next, the constructed 
mode was considered as an initial condition for Eq. 
and was allowed to evolve freely in time. The evolution 
is shown in Fig. O (c). 

In Fig. 131 (c) one can observe relatively stable dynam- 
ics, although accompanied by emission of some radiation, 
what means that the analytic approximation is relatively 
good. Meantime, any significant deviation of the poten- 
tial amplitude from the unity, has led to unstable dynam- 
ics of the initial conditions obtained from H37a|) . I)37b|) 
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FIG. 3: (a) Stationary coupled mode tp{x,0) for a = 0, P — 1 
and the corresponding envelopes ipo and ipi are shown, (b) 
Combined potential V{x) — Vt{x) + V-{x). The parameters 
are: A = -1, rj = 1.5, £ = 0.1, xa = 0.5, E = 0.9688 (eo = 
—El — 0.4982). Some parameters such as Waa, Ma are taken 
from Table I. (c) Time evolution of the coupled mode shown 
in (a). The solid black lines indicate the initial (top) and final 
(bottom) profiles. 



and (|24() . what is natural in view of the limitations of 
the theory described above. 



IV. DRIVING DEFECT MODES THROUGH 
THE OPTICAL LATTICE 

A. Response of defect modes to a simple shift of 
the defect position 

^From the physical point of view and also thinking 
on possible applications of defect modes, it is relevant 
to understand if they can be driven through the lattice 
and what are the limitations of those motions. Following 
the analysis of Ref. |0, here we start by studying the 
behavior of defect modes under a shift of the defect po- 
sition. We first consider evolution of the simplest cosh- 
and sinh- modes affected by the shift of the defect center 
by a lattice period (i.e. by tt in our case). As one can see 
from Fig. 0(a), (b) (for cosh-mode) and Fig. 21(e) i (f) 
(for sinh-mode), in this situation, the defect modes fol- 
low the position of the center of the defect potential, i.e. 
after some relaxation time the defect mode follows the 
defect position. However, when the defect is shifted by 
a half-period, both modes become unstable and spread 
out as it is shown in Fig. (c) , (d) (for cosh-mode) and 
Fig. (g), (h) (for sinh-mode). The observed behav- 
ior can be understood by realizing that a positive defect 
corresponds to local repulsion of the atoms having real, 
and thus positive, mass. Then, being placed on top of 
the maximum of the lattice potential the defect results 
in strong local repulsion of the atoms, what leads to de- 
cay of the defect mode. This does not happen when the 
shift of the defect mode is tt since in that case the local 
repulsion induced by the defect is compensated by the lo- 
cal attraction of the lattice potential, and for the whole 
excitation the effective mass approximation is applied. 

In the case of a two-hump cosh-mode the picture 
changes drastically. By shifting the defect to the next 
adjacent minimum (in our case to the right) we observe 
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FIG. 4: Dynamics of the defect modes under defect shifts in 
the positive direction at time t = 10 by tt (left column) and 
by 7r/2 (right column). The blocks of four panels from the 
top correspond to one- hump cosh-mode (a)-(d); sinh-mode 
(e)-(h); two-hump cosh-mode (i)-(l); one-hole tanh-mode (m)- 
(p); and two-hole tanh-mode (q)-(t). Initial profiles (dashed 
lines) and defect parameters in (a,c), (e,g), (i,k), (m,o) and 
(q,s) are the same as those in Fig. |21 (a), (c), (e), (g) and 
(i), correspondingly. The final profiles (solid lines) are taken 
from the final time of the corresponding density plots. In the 
panels with density plots the horizontal dashed line shows the 
time of the shift of the defect. 



that the mode breaks into two parts: a one-hump mode, 
which after a relaxation time stabilizes around the new 
potential minimum where the defect is centered, and a 
mode moving outwards the defect with a constant veloc- 
ity [see Fig. If; however, the defect is shifted 
to the nearest potential maximum (i.e. to Xd — 7r/2) 
the two-hump cosh-mode evolves into a one-hump mode 



stabihzed at the position of the defect [Fig. 0] (k) , (1)] . 
While this behavior may seem contradictory to what we 
observed in the case of the one-hump cosh-mode, it can 
be understood by taking into account that now the defect 
is negative: Vq < 0, and thus acts as repulsive impurity 
for the wavepackets having negative effective mass and 
as attractive impurity for real atoms. 

Let us now think of the two-hump solutions as two 
wavepackets whose attraction is balanced by the repul- 
sive defect placed between them. Then, a shift of the de- 
fect to the right by one lattice period destabilizes this bal- 
ance and pushes the right-placed wavepacket in the posi- 
tive direction, what leads to its motion [see Fig.^(i),(j)]. 
Also, the attraction between the two wavepackets results 
in a small shift of the left-placed wavepacket toward the 
defect. When the defect is shifted by 7r/2 and is placed 
on top of the nearest local maximum of the OL [see Fig. 
0](k),(l)], its effect is dramatically reduced by the lattice 
maximum, which acts in the opposite sense. Therefore, 
the repulsion between the left and right waves is reduced 
what leads to their motion toward each other and col- 
lapsing in one localized wavepacket. The proposed inter- 
pretation of the phenomenon was supported by further 
studies of the dynamics of two- hump cosh- modes. First, 
after a sudden switch off of the defect without any dis- 
location, we observed a merging of the two wavepackets 
very similar to one shown in Fig.^(l). Next, we have ob- 
served continuing motion of the wavepacket (after separa- 
tion with the defect, it can be more appropriately called 
envelope soliton) as it would be reflected by the defect, 
although weakened by the lattice potential. The first 
stages of such dynamics are shown in Fig.0](l). Finally, 
by increasing the strength of the defect we succeeded to 
create a more localized wavepacket moving faster. 

The response of the tanh-modes, given by the formula 
(|33|l . under a shift of the defect is shown in Fig0](m)-(t). 
The one-hole tanh-mode follows the defect shift to the 
nearest potential minimum accompanied by excitation 
of oscillating dynamics of the mode amplitude. When 
the defect is shifted to the nearest maxima the one-hole 
tanh-mode breaks into two modes which move outwards 
in opposite directions [see Fig. ^(o), (p)]. In order to 
understand the dynamics of a dark mode, in the same 
spirit as this was done above for the bright modes, we 
interpret the mode as a dark soliton interacting with 
a delta-impurity and recall the results on such interac- 
tion [221 ■ In the case at hand Vo < and Mq < 0, 
what in terms of the equation H28|) means that the im- 
purity attracts the soliton having initially zero velocity: 
as it is shown in Ref. J,2] [see the formula (16) there] the 
initial acceleration of the soliton due to attraction is pro- 
portional to MqVq. This explains the dynamics shown in 
Fig.2|(m), (n). When the defect has a tt/2 shift, as it was 
explained above, its strength in the final position is signif- 
icantly weakened (it is now placed on the top of a lattice 
maximum) . Since the shift of the defect acts as a pertur- 
bation of the initials conditions which leads to the decay 
of a static soliton into two solitons with finite velocities 
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[moving outwards the defect according to the momentum 
conservation law, as it is shown in Fig. 2|(o), (p)]. 

The decay of defect modes into two solitons is more 
evident in the case of a two-hole initial condition and 
with a repulsive impurity which is repulsive for the static 
soliton (see Eq. (16) in Ref. ^) shown in Fig|3|(s), (t) 
for the tt/2 shift of the defect position. At the same time, 
the two-hole tanh-mode displays more sophisticated evo- 
lution, shown in the panel (q), (r) of Fig. 0] After the 
shift of the defect to the nearest lattice minimum, the two 
holes start to behave as two dark solitons, one of them 
(the one closer to the new position of the impurity) being 
trapped while the other one moves far from the impurity. 



B. Driving the defect modes 

In Ref. 10] it was shown that the robustness of the 
cosh-mode described in the preceding subsection allows 
one to drive the defect mode along the lattice. Since the 
defect mode is destroyed when it stays close to the local 
maxima of the potential, in order to provide stable mo- 
tion of the mode, the defect motion cannot be uniform, 
but described by a step-like function which consists of 
two characteristic time intervals: a fast enough one, t, 
in which the defect is shifted by one lattice period, and 
another larger one, T, allowing for the defect-mode to 
recover its shape on the next lattice site. 

The physical reasons for the existence of these two 
characteristic times, r and T are rather simple. Indeed, 
the respective microscopic dynamics can be interpreted 
as a tunneling of part of atoms, attracted by the shifted 
defect, from the potential well x = to the neighbor 
potential well x = tt. Let the characteristic tunneling 
time be tq. Taking into account that all intermediate 
states (i.e. those with center different from 7rn, where 
n is an integer) are unstable, one finds that the require- 
ment r <C To <C T must be satisfied. This last constrain 
evidently limits the speed with which the mode can be 
driven along the lattice. 

In order to find a numerical estimate for the tunnel- 
ing time To we take into account that we are working 
in the framework of the effective mass approximation 
which implies relatively small differences in populations 
of neighbor potential wells. Also we are dealing with the 
modes bordering the boundary of the BZ, what means 
that the phase difference between atoms in adjacent cells 
is TT. Thus we can use the results for the half-period 
of oscillations of the linear tunneling of condensate in a 
double- well potential [23l |: 



To 



with 



is: = - 



dwQo dwpi 
dx dx 



AwqqWqi cos(2x) 



(40) 



dx, (41) 



for a particular potential given by H27() (here we have 
taken into account that the WFs make up an orthonormal 
basis). In particular, for A ~ —\, the mobility conditions 
are that the defect must be shifted with a velocity much 
bigger than 



3.67 



0.86, 



and then leave the system to relax for a time 
T > To « 3.67. 



(42) 



(43) 



In Ref. (see Fig. 4) we have presented an example 
of how the defect cosh-mode can be driven along the lat- 
tice by moving the position of the defect. According to 
the properties shown Fig. ^(e), (f), the sinh-mode also 
can be moved through the lattice. An example is shown 
in Fig|31 where we present the dynamics of the mode as 
well as the trajectory of the center of the wavepacket (x) 



(x) = ^ j x\i;\^dx 



and of its dispersion 



(44) 



(45) 



Here as usually N — J j^pdx is a number of particles. 

Due to the higher scattering on the potential barri- 
ers after each defect, shifting the dynamics of this mode 
is less robust than for the cosh-mode what leads to its 
faster dissipation as shown in Fig|5| Also, due to sig- 
nificant reflection on each period of reconstruction one 
observes larger deviation of the center of mass of the 
mode (x) (t) from the defect trajectory Xd{t) as well as 
significant growth of the width of the wavepacket D [see 
FigEl(a), (b)]. Therefore in the following analysis of the 
motion of defect modes through the lattice we will use 
cosh-mode. 

If the motion of the defect satisfies the conditions H42|l 
and (|43|l . more sophisticated dynamics can be generated 
by moving defects. For instance, in Fig. we show a 
zigzag trajectory of the cosh-mode induced by the corre- 
sponding zigzag motion of the defect. The defect mode 
experiences some changes when driven through the lat- 
tice. For instance, when driven over ten lattice periods 
the amplitudes of the initial and final maxima differ by 
about 10%. 

We have also carried out numerical experiments on col- 
lisions of defect modes driven by counter-propagating de- 
fects. In Fig. [7|we show an example with two cosh-modes 
each of them shifted by ten periods from the center in 
opposite directions. By moving the defects towards each 
other we induce a head-on collision at t = 500. When 
the defects continue to be driven after the collision [Fig. 
[3a)-(c)] they pass through each other and emerge with 
substantial shape modifications but still being identifi- 
able [Fig. EJc)]. When the defects are stopped at the 
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FIG. 5: Density plot of the sinh-mode having the initial pa- 
rameters as in Fig. |21(c) driven by the moving defect following 
trajectory Xd{t) = T^J2j^o ^(^ ~ j) ^^^^ ^ = 50. The solid 
black lines indicate the initial (top) and final (bottom) pro- 
files. In the inset (a) the dynamics of the center of mass (x) 
of the wavepacket (thick black line) and of the defect position 
Xd(t) (thin blue line) are shown. In inset (b) the dynamics of 
the wavepacket dispersion D(t) is presented. 
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FIG. 6: Density plot of driving of a cosh-mode with initial 
parameters as in Fig. |5| (a) by moving the defect following a 
zigzag trajectory according to the law Xd(t) — ttX] =o ~ 
to — T j) along each direct interval with T — 50. At times 
t = 250 and 650 we invert the direction of the motion of the 
defect which is followed by the inversion of the direction of 
the motion of the mode. The solid black lines indicate the 
initial (top) and final (bottom) profiles. In (a) dynamics of 
the center of mass {x){t) of the wavepacket (thick black) is 
shown. The defect position Xd{t) is indicated by the thin blue 
line in the same plot. Inset (b) presents the dynamics of the 
wavepacket width D{t). 



collision time t ~ 500 a single defect mode with larger 
amplitude and stronger localization appears [Fig. e), 
(f)]- 
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FIG. 7: Two regimes of collision of two identical cosh-modes, 
each having the initial parameters as in Fig. |2| (a). The 
initial distance between the modes is 20 periods. In both 
columns initially the modes are driven by the defects toward 
each other according to the law used in Fig. |S| and |S| with 
T = 50 what results in a collision ai t = 500 (the time where 
the defect positions coincide) . In the left column we show the 
dynamics where after collision defects continue to drive modes 
along the initial directions. In the right column the results 
correspond to the situation when the defects are stopped after 
they are met. Panels (b) and (e) show the respective temporal 
dynamics, while panels (c) and (f) show the final outputs at 
t = 10^ 



V. GENERATION OF DEFECT MODES 

A relevant question for experimental studies of defect 
modes is their generation in real situations. One of the 
possibilities consists in generating, first, a gap soliton 
with energy near a gap edge in an homogeneous OL 
and subsequently adiabatically switching on the defect by 
changing the defect amplitude from zero to some value. 
This will lead to the transformation of a gap soliton into a 
defect mode. A numerical implementation of this scheme 
is illustrated in Figs. ISjandO where we have presented 
the dependence of the number of particles N vs the en- 
ergy of the mode E for the original gap soliton (thin solid 
lines) and the emerging defect modes (thick dashed and 
dotted lines), found by searching stationary modes of Eq. 

More specifically, in order to generate cosh- and two- 
hump cosh-modes we start with the envelope soliton (as 
an example consider, say, the point I in Fig|SJ) with en- 
ergy E = Eq + Eq. Next, by increasing rj [the defect shape 
is chosen to be ()35|l ] we generate one-hump cosh-mode, 
shown as I^C transition [in Fig. |Hl(b) and|H|(b)], and by 
decreasing rj (i.e. increasing \r]\) we generate a two- hump 
cosh-mode, shown as (I^B) transition. In both cases we 
used the stationary ID GP equation with a periodic po- 
tential to compute numerically the stationary solutions 
and then solving the time-dependent equation JQ) with 
adiabatic growth of the defect amplitude according to 
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FIG. 8: [Color online] Adiabatic excitation of cosh- and two- 
hump cosh-modes is shown for A = —1, a = 1 and I = 0.1. 
Different branches N{E) in the first gap correspond to the 
lowest gap soliton for rj = (black, solid line), cosh- mode 
for 77 = 0.1 (red, dashed line) and two-hump cosh-mode for 
rj = —0.1 (blue, dotted line). In (a) and (b) the initial profile 
of the gap soliton (black, solid line) corresponds to the point 
I and final profiles of the two-hump cosh- and cosh-modes 
(blue dotted and red dashed, thick lines) in the points B and 
C, correspondingly, are shown which coincide with the profiles 
of the defect modes obtained by adiabatic growth of the defect 
amplitude. 

the law 77 — ±0.1 8111^(10^"^ 7rt). It worth pointing out 
that in the described transitions the numbers of particles 
was preserved. 

We have followed the same procedure for the genera- 
tion of a sinh-mode as schematically shown in Fig. |51 
However, the difference from the previous case is that 
to generate a sinh-mode from the gap soliton not only 
the defect amplitude must be increased adiabatically but 
also the sign of the interparticle interactions, a, must be 
changed from repulsive to attractive. Due to the fact 
that gap soliton does not exists for attractive interaction 
it starts to spread and the only small number of parti- 
cles transforms into the sinh-mode [see the process I^B 
in Fig. The resulting number of particles and corre- 
sponding profile of the sinh-mode slightly depends on the 
velocity of growth of the defect amplitude [in the present 
example we used 77 = 0.5sin^(5 x 10~^7rt)]. 



VI. CONCLUSIONS 

In this paper we have investigated localized defect 
modes of the one-dimensional Gross-Pitaevskii equation 
in the presence of an optical lattice and a localized defect. 

Using the expansion of the condensate wave function 
over Wannier functions it was shown that the problem 
can be reduced to the study of the dynamics of the en- 
velope which is governed by the nonlinear Schrodinger 
equation with a delta-impurity (in the one-band approx- 



FIG. 9: [Color online] Adiabatic excitation of a sinh- and 
cosh-mode mode is shown for A = —1, r] = 0.5 and £ = 0.1. 
Different branches N{E) in the first gap correspond to the 
lowest gap soliton for = and a — 1 (black, thin line) , cosh- 
mode for ri = 0.5 and a = 1 (red, dashed line) and sinh-mode 
for -q = 0.5 and a = —1 (blue, dotted line). In (a) and (b) the 
initial profile of the gap soliton (black, solid line) corresponds 
to the point I and final profiles of the sinh- and cosh-modes 
(blue, dotted and red, dashed thick lines) in the points B and 
C, correspondingly, are shown which coincide with profiles of 
the defect modes obtained by adiabatic growth of the defect 
amplitude. 

imation), or in a general case by a system of coupled 
nonlinear Schrodinger equations with delta-impurities (in 
the two-band approximation). The symmetric analytical 
solutions the obtained equation in the stationary case de- 
scribe envelopes of the stationary defect modes of cosh — , 
sinh—, and tanh— types. Our approximate analytical 
formulas for impurity modes are in good agreement with 
results obtained by direct numerical calculations. 

We have also verified the stability of the impurity 
modes by direct numerical simulations of the underlying 
one-dimensional Gross-Pitaevskii equation with a peri- 
odic potential and a localized defect showing that the 
defect modes are indeed stable under small perturbations 
(i.e. linearly stable). 

Using the same ideas developed in Ref. we have 
shown that impurity cosh-modes can be made to follow 
complex trajectories of the defect position while sinh- 
modes cannot be driven so easily. An example is shown 
where a zigzag motion of a cosh-mode is induced. We 
have also used this control to induce collisions of two 
mutually displaced defect modes. The outcome depends 
on the defect trajectory after the collision being either a 
pair of driven defect modes or a single defect mode with 
stronger localization. 

Finally we have discussed in detail how to excite differ- 
ent types of these defect modes in realistic scenarios. We 
think that the mechanism proposed in this paper could 
provide one more possibility for experimentalists to man- 
age gap solitons by using a small localized defect of an 
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optical lattice. 
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